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1.0  INTRODUCTION 


In  a  band  of  rather  low  frequencies,  almost  all  of  the 
ambient  noise  in  the  ocean  is  generated  by  ship  traffic.  The 
performance  of  sonars  operating  in  this  band  is  limited  by 
this  component  of  background  noise;  therefore,  it  is  necessary 
to  measure  or  predict  certain  statistical  measures  of  beam  noise 
to  support  system  design  and  deployment  decisions.  ^ 

^  The  probability  measures  of  beam  noise  depend  on  many 
factors:  array  configuration,  orientation,  and  location 
(including  depth),  and  season.  Thus  it  would  be  both  expensive 
and  time-consuming  to  attempt  to  obtain  the  required  beam  noise 
measures  for  a  large  number  of  sites  experimentally.  Further¬ 
more,  if  shipping  routes  and  traffic  densities.  wen5"to  change 
because  of  strategic  or-^oun'amiC'^cSuses,  the  measurements  would 
have  to  be  repeated. jlThere  is,  therefore,  considerable  incentive 
for  the  development  of  a  capability  of  predicting  statistical 
measures  of  array  beam  noise. 


^^>The  aspects  of  beam  noise  that  should  be  considered  depend 
on  the  type  of  performance  predictions  desired.  For  detection 
predictions,  the  quantity  of  interest  is  the  total  noise  power 
in  selected  frequency  bands.  However,  for  prediction  of  false 
alarm  or  false  classification  rates,  characterization  of  the 
narrowband  components  of  shipping  noise  must  be  considered. 

Both  aspects  are  addressee^  in  this  report y  however,  the  bulk 
of  the  report  pertains  to  the  former. 
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The  development  of  a  beam  noise  prediction  capability  is 
partly  based  on  prior  work  reported  in  Reference  1,  which 
concerned  prediction  of  noise  at  a  receiving  point  due  to 
ships  in  a  sector  of  azimuth.  Reference  2  introduced  concepts 
and  approaches  that  are  more  completely  covered  in  Reference  1. 

Certain  input  data  required  for  the  predictions  are  given 
in  other  reports.  Characteristics  of  noise  radiated  by  surface 
ships  are  given  in  References  3>  4,  6  and  7-  Data  on  ocean 
route  envelopes  are  given  in  Reference  5* 

Sections  2  through  4  concern  the  prediction  of  the  prob¬ 
ability  density  of  the  random  variable  Y,  the  total  averaged 
noise  power  in  a  specified  passband  at  the  output  of  a  beam- 
former.  The  approach  is  to  first  calculate  the  characteristic 
function  of  Y,  and  then  to  transform  it  to  the  density  function. 

Section  2  presents  an  expression  for  Y,  and  derives  a 
formula  for  the  calculation  of  the  characteristic  function. 

It  also  defines  a  geometry  for  segmented  routes,  and  derives 
a  function  that  characterizes  the  distribution  of  ship  coordi¬ 
nates  normal  to  those  routes. 

Section  3  concerns  the  calculation  of  the  probability 
density  of  Y  from  its  characteristic  function.  Since  Y  is 
non-negative,  it  is  shown  that  its  density  can  be  calculated 
from  one  of  two  alternative  formulas,  each  involving  only  a 
single  integral.  A  special  numerical  transformation  technique 
is  then  derived  for  one  of  the  alternative  formulas. 
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Section  addresses  the  input  data  for  the  calculation  of 
the  characteristic  function  of  Y,  including  a  transmission  loss 
function,  a  beam  pattern  function,  and  the  geometry  linking  a 
specific  sensor  to  the  shipping  routes  in  the  surrounding  ocean. 
A  means  is  also  devised  and  utilized  for  obtaining  source  charac 
teristic  functions  for  various  types  of  ships.  These  functions 
are  required  for  the  calculation  of  the  characteristic  function 
of  Y. 


Section  5  derives  a  formula  for  calculating  the  expected 
value  of  N,  a  random  variable  representing  the  number  of  lines 
observed  at  a  beamformer  output  whose  center  frequencies  fall 
within  a  specified  frequency  interval,  and  whose  average  power 
exceeds  a  specified  threshold. 


A  companion  report.  Bolt  Beranek  and  Newman  Inc.  Report  No. 
3654,  is  entitled  "A  Computer  Program  for  Beam  Noise  Prediction. 
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2.0  CHARACTERISTIC  FUNCTION  FOR  Y 

This  section  derives  a  formula  for  the  characteristic 
function  for  the  total  noise  power  in  a  specified 

band  at  the  beamformer  output.  It  also  deals  with  related 
matters  such  as  a  geometry  for  segmented  routes,  and  a  func¬ 
tion  that  characterizes  the  distribution  of  ship  coordinates 
normal  to  those  routes.  First  however,  the  random  variable 
Y  will  be  stated  in  terms  of  the  contributions  of  individual 
ships  exposed  to  the  sensor. 

2.1  Derivation  of 

The  averaged  noise  power  at  the  output  of  a  beamformer 
in  a  specified  passband  resulting  from  ship  traffic  in  an 
acoustic  basin*  may  be  expressed  as 


Y  = 


A., 

m  n  ij 
III 
i=l  J=1  k=l 


SijkZijkBijk 


(1) 


where  m  is  the  number  of  routes  in  the  basin. 


n  is  the  number  of  ship  types. 


is  the 
^  Aij 


number  of  ships  of  type  j 
is  a  random  variable, 


on  route 


*The  acoustic  basin  in  which  the  array  is  located  is  defined  by 
bathymetry  (e.g. ,  land  masses  and  high  underwater  ridges.)  The 
intensity  of  sound  from  sources  outside  of  the  basin  is  negli¬ 
gible  . 
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S. ..  is  the  source  intensity  of  the  k —  ship  of 

1 J  K 

type  j  on  route  i;  it  is  a  random  variable 
that  is  statistically  independent  of  the 
source  intensity  of  any  other  ship. 


7 

ijk  is  the  intensity  transmission  ratio  from 
ship  ijk  to  the  receiving  point, 

Bij<k  is  the  gain  for  a  plane-wave  arriving  at 
the  array  from  ship  ijk. 


Section  1.2  of  Reference  1  shows  that  (1)  is  a  consequence  of 
the  principle  of  superposition  of  the  instantaneous  pressures 
of  the  sound  waves  from  multiple  point  sources. 


The  probability  density  function  for  Y  can  be  obtained 
from  its  characteristic  function 

♦y  (w)  =  E[expJ(i)Y]  (2) 


where  E  denotes  the  expectation  operator.  Substituting  (1) 
in  (2)  gives 


4>y(w) 


SlJkZiJkBiJk 


Sm  n  AiJ 

expjw  £  £  £ 

i=l  J=1  k=l 

|*  m  n  AiJ  1 

E  it  IT  IT  expjwSijKZ  kBiJk 

[i=l  J-l  k-1  J 


IT  IT  4>  .M 
1=1  J=1  10 


(3) 

(4) 

(5) 
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where 


riJ 


(w) 


AiJ 

k!j. 


(6) 


In  the  sequel  it  is  assumed  that 


Zijk  =  z(Gijk*  Qijk) 
Bijk  =  t(Gijk*  Qijk5 


(7) 

(8) 


where  z(  ,  )  is  a  deterministic  transmission  ratio  function; 
b(  ,  )  is  a  deterministic  gain  function  for  a  plane-wave;  GTT„ 
is  the  earth-centered  angle  measuring  the  longitudinal  coordi¬ 
nate  (in  the  direction  of  the  nominal  route)  of  ship  ijk;  and 
Qijk  the  earth-centered  angle  measuring  the  transverse  coor¬ 
dinate  (normal  to  the  nomial  route)  of  ship  ijk.  It  is  assumed 
that  all  ship  coordinates  are  random  variables  that  are  mutually 
independent.  Substituting  (7)  and  (8)  in  (6)  and  performing 
the  expected-value  operation  gives 


*  (w)  =  r  p  (a)[/dgf  (f,)/dqf  (q;g)/dsf  (s) 

1J  a=0  lJ  L  Gi  Qi  SJ 

5,q)b(G.q)J 


x  expju«z(g. 


(9) 
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where  p^ta)  =  p(Ajj  =  a  =  0,  1,  2,  3,  ... 

fQ  (  )  is  the  probability  density  function  for 
Gijk*  J  =  15  2>  k  =  1»  2’  3’  "• 


fn  (  ,  )  is  the  probability  density  function  for 
^i 

Qijk  dePendinS  on  the  variable  g;  j  =  1,  2, 
•  •  •  ,  n j  k  1 ,  2 ^  3,  ... 


The  inner  integral  of  (9)  can  be  expressed  as  a  characteristic 
function;  thus 

00 

♦ijlu)  =  P^U)  I  fdgFGi(E)/dqfQi(q;g) 

x  4>gJ  [to  z(e.q)b(g,q)))a  (10) 

where  (  )  is  the  source  characteristic  function  for  ships  of 
type  j. 

In  the  sequel,  it  is  assumed  that  A ^  is  a  Poisson  random 
variable;  thus,  its  probabilities  are  given  by 

PjjU)  =  exp(-aij)aiaj  v  a!,  a  =0,  1,  2,  3,  ...  (H) 

where  a^  is  the  average  number  of  ships  of  type  j  in  an  in¬ 
crement  of  route  i.  Substituting  (11)  into  (10)  gives 
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♦jjlw)  =  expl-a^)  £  (aiJ/dgfGi(g)/dqfQi(q;g) 

x  4>gj  Iwz  (g,q)b(g,q)]}a  a.'  (12) 

The  sum  is  the  power  series  expansion  for  the  exponential  of 
the  quantity  in  braces;  thus 

*^("0  =  exp(-aiJ)  exp{aiJ/dgf(Ji(g)/dqfQi(q;g) 

x  4>Sjlwz  (g,q)b(g,q)]}  (13) 

The  cumulant  function  is  the  logarithm  of  the  characteristic 
function 

Vw)  = 

x  *SJ  z  (8»q)b(g,q)]  -1}  (1*0 

If  the  distributions  of  the  longitudinal  coordinates 
are  uniform  over  the  route,  then 

fGi(g)  =  t'*  0  1  «  £li 

=0,  elsewhere  (15) 

where  V,±  is  the  length  of  route  i.  Substituting  (15)  in  (lM) 
yields 
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f 

,(w)  =  pk  /i  dg/dqf  . (q,c) 

1J  1J  Q  H1 


4>Sjluz  (e,q)t>(G»q)l  -  p/ikij 


(16) 


where  p  Is  the  earth's  radius,  and  is  the  average  number  of 
ships  of  type  j  per  unit  length  of  route  i.  In  (16)  the  variables 
of  integration  are  the  route  variables  g  and  q;  an  alternative 
set  of  variables  is  the  pair  of  observation  variables  related  to 
the  range  and  bearing  of  a  point  from  the  sensor.  These  variables 
are  defined  in  Section  2.2,  and  the  alternative  to  (16)  is  given 
by  (33). 


The  characteristic  function  is  the  exponential  of  the  cumu- 
lant  function;  thus,  (w)  =  exp  1^(10).  Substituting  this 
result  in  (5)  gives  the  characteristic  function  of  Y 


m  n 

$v(w)  =  IT  IT  exp  41.  ,(w)  (17) 

1  i=l  J=1  J 

m  n 

=  exp  Z  Z  <('• .  («) 
i=l  J=1  J 


The  second  of  these  equivalent  forms  is  easier  to  calculate. 


2.2  Geometry  for  Beam  Noise  Prediction 

Section  2.2  of  Reference  1  discusses  a  coordinate  system 
that  is  applicable  to  the  prediction  of  characteristics  of  noise 
at  a  point  received  from  ships  in  a  narrow  sector  of  observa¬ 
tion. 
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The  prediction  of  beam  noise,  however,  requires  accounting 
for  sidelobe  contributions  as  well  as  the  main  lobe  contribution. 
As  a  consequence,  the  observation  sector  for  a  given  route  is 
determined  by  the  extent  of  the  route  in  the  basin  that  includes 
the  observation  point.  A  further  complication  is  that  within 
a  basin  a  route  may  consist  of  several  great  circle  segments. 

For  these  reasons  it  is  important  to  devise  a  coordinate  system 
that  is  suitable  for  the  efficient  calculation  of  beam  noise. 

Figure  1  illustrates  a  route  with  three  segments  with 
end-points  through  R^.  The  observation  point  is  designated 
by  R  .  This  case  defines  the  three  spherical  triangles  shown 
in  Figure  1,  which  also  defines  arcs  and  angles  that  will  be 
required. 

The  calculation  of  will  be  executed  in  spherical 

polar  coordinates.  One  of  the  coordinates  is  the  range  arc  r, 
the  earth-centered  angle  between  the  observation  point  and  a 
ship.  For  each  route  segment  i,  an  angle  is  defined  which 
is  the  angle  between  the  arcs  si  and  r.  These  angles  will  be 
discussed  in  greater  detail. 

The  calculation  of  the  parameters  required  for  the 
calculation  of  $>y(u))  can  be  organized  by  means  of  a  data 
array  whose  headings  are  shown  in  Table  1,  along  with  explana¬ 
tory  notes  and  the  numbers  of  figures  in  which  the  parameters 
are  defined.  The  order  of  data  generation  is  from  top  to  bottom. 
Item  0  is  the  point  number,  that  is,  the  subscript  of  R^.  As 
shown  in  Figure  1,  Rq  is  the  observation  point,  and  points  R^ 

R 2,  R^,  ...,  Rn+1  are  the  numbers  of  end-points  of  the  n  route 
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Table 

:  1:  Parameters  for  a  Route 

Item 

Heading 

Comments 

0 

i 

subscripts  of  Ri,  i  =  0,  1,  2,  3, 

n.  Figure  1 

1 

latitude 

of  point  Ri,  input  data 

2 

longitude 

of  point  Ri ,  input  data 

3 

ci 

(1) ,  Figure  2 

4 

Ni 

(2) ,  Figure  2 

5 

COS  Si 

(3) ,  Figure  2 

6 

Ci 

(4) ,  Figure  2 

7 

zi 

azimuth  of  Ri,  Eq.  (5) 

8 

(6) ,  Figure  1 

9 

sin  li 

(7) ,  Figure  1 

10 

Fi 

(8) ,  Figure  1 

t 
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segments.  The  convention  for  numbering  the  end-points  is  as 
follows:  from  a  point  within  the  spherical  polygon  Rq,  F^, 

R^j  .  ...,  ^n+i >  Rj.  is  the  route  end-point  that  is  the  more 
counter-clockwise  around  Rq. 

Items  1  and  2  of  Table  1  are  input  data,  the  latitudes 
and  longitudes  of  the  points.  Items  3  and  defined  in 
Figure  2,  are  calculated  from  Items  1  and  2.  The  sides  of  the 
(north)  polar  triangles  are 

Ci  =  tt/2  -  north  latitude 

=  tt/2  +  south  latitude  (18) 


The  polar  angles  are  given  by 

Ni  =  lWi  "  wol  (19) 

where  Wi  =  west  longitude 

=  2tt  -  east  longitude 

Items  5  and  6  are  obtained  by  solving  the  polar  triangles 
cos  =  cos  c^  cos  Cq  +  sin  c^  sin  cQ  cos  (20) 

sin  C.  =  sin  c.  sin  N.  sin  s.  (21) 

l  li  l 
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Item  7  Is  derived  from  Item  6.  The  azimuth  of  point  is 


z . 
1 


c. ,  W.  <  Vtf 

l  1—0 


2tt  -  C.  ,  W.  >  W. 

l  i  0 


(22) 


where  VT  is  given  below  (19).  The  internal  angles  (Item  8)  of 
the  route  segment  triangles  shown  in  Figure  1  are 


(23) 


Items  9  and  10  are  obtained  from  solutions  of  the  route 
segment  triangles 


cos  1^  =  cos  s^  cos  s^+1  +  sin  s^  sin  s^+1  cos  1^  (24) 

sin  Fi  =  sin  si+1  sin  Ii  4-  sin  li  (25) 

It  might  be  useful  to  add  the  item  sin  si  to  Table  1 
since  this  quantity  is  required  for  several  calculations. 


The  observation  variables  in  route  segment  i  are  shown  in 
Figure  3.  The  variable  ri  is  the  arc  from  the  observation 
point  to  the  ship,  and  is  the  surface  angle  between  the 
arcs  s^  and  ri-  The  ship  position  is  also  given  by  the 
variables  g.^  and  q^ ;  the  former  is  the  arc  F^T  in  the  direction 
of  travel,  and  the  latter  is  the  transverse  arc  TS.  The 
relationship  between  these  pairs  of  variables  is  required. 
Auxiliary  variables  involved  in  the  derivation  are  obtained 
from  solutions  of  the  triangle  Rq,  S: 
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cos  cL  =  cos  s^  cos  +  sin  s^  sin  cos  (26) 

sin  R^  =  sin  sin  *  sin  (27) 

where  Ri  is  the  surface  angle  between  arcs  s^  and  d^.  The 
variables  g^  and  qi  are  found  from  solutions  of  the  right 
triangle  FLST: 


tan  g.  =  cos 

(R.  -  F. ) 

tan  d . 

(28) 

l  l 

l 

sin  =  sin 

(R.  -  F.) 
l  1 

sin  d . 
i 

(29) 

For  each  route  segment,  the  following  numerical  inte¬ 
gration  regime  can  be  considered.  The  variable  Di  is  advanced 
in  increments  from  zero  to  1^,  and  for  each  selected  value 
of  D^,  the  variable  r^  is  stepped  through  the  appropriate 
range,  determined  per  the  following.  For  a  ship  that  is  on 
the  nominal  track,  the  points  S  and  T  in  Figure  3  coincide, 
as  do  d^  and  g^.  The  first  value  of  r^  used  is  found  from 
the  simultaneous  solution  of 

tan  ^ (r\  -  =  tan  h  s ^  sin  -  D^)  t  sin  (F±  +  Di) 

(30) 

tan  h  (r .  +  g . )  =  tan  h  s.  cos  *s(F.  -  D.)  t  cos  %(F.  +  D.) 

(3D 
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The  variable  r^  is  increased  by  increments  until 

IgJ  1  h  wi(gi}  (32) 

where  w^g)  is  the  width  of  the  route  at  Then  starting 

again  at  the  initial  value  of  r^ ,  the  value  is  decreased  until 
(32)  is  again  satisfied. 

The  procedure  described  does  not  completely  characterize 
the  routes  at  the  turns,  and  partially  self-cancelling  errors 
will  occur  near  the  segment  end-points.  This  problem  could  be 
remedied  in  a  straightforward  fashion  with  the  expense  of 
greater  complexity,  which  does  not  appear  warranted  with  the 
current  state  of  the  art. 

In  the  previous  section,  a  formula  used  in  the  calculation 
of  the  characteristic  function  of  Y  employed  the  ship  route 
variables  in  the  integrals.  This  formula  may  now  be  restated 
in  terms  of  the  observation  variables  r  and  D. *  The  alternative 
to  (16)  of  Section  2.1  is. 

(w)  =  pkij  ^  dr  !  dD|j(r,D)|  fq  (g,q)  S>Sj  [wzt  (r ,  D)  (r ,  D)  ] 
-P4kij  (33) 

where  the  Jacobian  is  given  by  (2-23)  of  Reference  1. 


#To  simplify  the  result,  the  route  segment  subscript  has  been 
suppressed. 
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A  recently  distributed  report  [51  presents  certain  data 
characterizing  ocean  route  envelopes:  the  latitude  and  longi¬ 
tude  of  the  end-points  of  route  segments;  the  (transverse) 
width  of  the  segments  at  their  end-points;  and  the  maximum 
width  of  the  segment.  The  width  of  a  route  may  be  given  as 
a  function  of  the  distance  travelled  along  a  segment.  A  quad¬ 
ratic  formulation  is 

Vgi}  =  wi  +  b.g.  +  e.g*  (34) 

where  bi  =  2  £,7* 

ei  “  ^  [""i 


w^  is  the  route  width  at 


nr  is  the  maximum  width  of  the  segment  R^R^+1 


The  probability  density  function  of  the  transverse  coordinates 
may  be  expressed  as  a  function  of  the  variable  and  the 
(variable)  parameter  wi(gi). 

The  next  section  derives  a  function  to  represent  the  distri¬ 
bution  of  ships'  coordinates  normal  to  the  nominal  route. 
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2.3  Ship  Distribution  Across  a  Route 


The  computation  of  the  characteristic  function  $y(w)  re¬ 
quires  a  probability  density  function  for  the  ship  position 
across  the  route.  In  Equation  (33)  of  Section  2.2  this  density 
function  is  denoted  for  the  i5il  route  by  fQ^Cq),  where  q  is  the 
earth-centered  angle  in  radians  representing  across-route  ship 
position  measured  from  the  center  of  the  route.  This  section 


develops  a  function  for 
density  function. 


(q)  based  on  the  beta  probability 


Define  the  width  of  the  route,  at  a  given  point  along  the 
route,  as  W,  where  W  is  an  earth-centered  angle  in  radians. 
Define  the  half-width  of  the  route  as 


qQ  =  (1/2) W 


(35) 


The  probability  of  a  ship  that  is  using  this  route  being  outside 

the  interval  -qQ  <  q  <  +  qQ  is  assumed  to  be  zero. 

Since  exact  information  of  ship  distributions  across  the 
world’s  shipping  routes  is  not  presently  available,  a  probability 
density  function  must  be  chosen  to  model  the  uncertainty  in  ship 
position  within  the  route  width  W.  An  appropriate  choice  is 
felt  to  be  the  beta  probability  density  function  with  the 
parameter  values  both  set  equal  to  two.  This  density  function 
is  symmetric  about  its  mean  value,  which  if  chosen  as  the  center 

of  the  route,  does  not  favor  one  side  of  a  route  over  the  other. 


i 
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The  density  function  to  be  used  for  across  the  route  ship 
position  uncertainty  is 


Where  qQ  is  as  defined  in  Equation  ( 35 )  ♦  Figure  4  shows  a  plot 
of  the  density  function  of  Equation  (36). 
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FIGURE  4  Density  of  Position  Across  Route 
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3.0  PROBABILITY  DENSITY  FUNCTION  FOR  Y 

The  probability  density  function  of  Y  is  calculated  from 
its  characteristic  function  4>y(oj),  which  in  turn  is  obtained 
using  (33)  of  Section  2.2  in  (17)  of  Section  2.1.  Since  the 
calculation  is  to  be  done  by  numerical  methods,  consideration 
must  be  given  to  selecting  a  set  of  values  of  co  at  which  is 
to  be  calculated,  and  to  an  efficient  method  to  achieve  its 
transformation . 

Section  3-1  considers  several  factors  pertaining  to  a 
uniform  sampling  scheme.  The  first  factor  considered  is  the 
reduction  of  the  number  of  sample  values  that  have  to  be  cal¬ 
culated.  It  is  shown  that  a  technique  similar  to  basebanding 
in  communications  systems  is  not  effective  if  the  random 
variable  is  non-negative.  For  such  variables,  however,  it  is 
shown  that  alternative  forms  are  available  for  calculating  the 
density  function  from  the  characteristic  function.  Each  of 
these  forms  involve  only  a  single  integral  vice  the  two  required 
for  a  variable  not  so  restricted.  Another  way  to  reduce  the 
number  of  samples  required  is  to  employ  higher-order  function 
interpolation  schemes  for  the  density  calculation.  A  second- 
order  scheme  is  selected  to  provide  a  reasonable  compromise 
between  sample  number  reduction  and  calculation  complexity. 

The  final  factor  considered  is  a  means  of  stabilizing  the  range 
of  integration  in  the  calculation  of  the  density  function. 

Section  3.2  develops  an  algorithm  for  calculating  the 
density  function  of  Y  from  the  real  part  of  the  characteristic 
function.  The  method  developed  is  attractive  in  that  the  selec¬ 
tion  of  the  inter-sample  interval  does  not  depend  on  the  value 
of  the  probability  density  variable  y.  This  would  not  have 
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been  the  case  had  a  conventional  numerical  integration  scheme 
been  employed. 

3.1  Function  Sampling 

The  previous  version  of  the  ambient  noise  prediction  program 
calculated  the  characteristic  function  of  Y  at  frequencies  uk 
that  were  exponentially  spaced.  There  are  several  disadvantages 
that  are  inherent  in  this  sampling  procedure. 

•  No  calculation  is  made  for  w  =  0  (This  would  be 
useful  as  a  check  point.) 

•  The  space  between  u>  =  0  and  the  smallest  value 
of  to  is  much  larger  than  the  space  between  it 
and  the  next  larger  value. 

•  The  space  between  adjacent  values  of  oj  at  the 
upper  range  is  so  large  that  the  approximation 
of  the  function  as  a  straight  line  between  the 
points  lead  to  large  errors. 

This  section  concerns  the  investigation  of  an  alternative  sampl¬ 
ing  regime. 

One  factor  to  consider  is  the  reduction  of  the  required 
number  of  samples  by  a  method  that  is  analogous  to  the  practice 
of  heterodyning  a  bandpass  signal  to  baseband.  The  general 
possibilities  are  illustrated  by  means  of  two  examples. 

Consider  the  possibility  of  calculating  the  characteristic 
function  for  a  standardized  random  variable 

S  =  (X  -  m)  f  o 


2k 


(37) 
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where  x  is  any  random  variable 

m  is  its  mean  value 

o2  is  its  variance 

It  is  easily  shown  that  the  characteristic  function  for 
the  standardized  random  variable  is 

4>s(to)  =  exp (- jwm  *  a)  $x(u  *  o)  (38) 

where  4>x(w)  is  the  characteristic  function  for  X. 

Suppose  X  is  Gaussian.  The  characteristic  functions  for 
it  and  the  standardized  variable  are  respectively 

4>x(w)  =  (expjmw)  exp  -  is  (aw)  2  (39) 

4>s  (oj)  =  exp  -  igw2  (40) 

It  is  apparent  from  (39)  that  the  modulation  frequency  of  the 
first  factor  is  proportional  to  the  mean  value;  the  number  of 
samples  required  for  accurate  sampling  is  thereby  increased. 

The  spread  of  the  second  factor  varies  inversely  with  the 
variance.  Thus  the  sampling  problem  induced  by  the  first 
factor  increases  with  a  decrease  in  variance.  If  the  mean 
value  of  X  were  zero,  it  is  surmised  that  the  number  of 
samples  required  is  independent  of  the  variance  if  the  sample 
spacing  is  proportional  to  the  spread.  Prom  (40)  it  is  seen 
that  the  characteristic  function  for  the  standardized  variable 
is  independent  of  the  distribution  parameters  of  X,  i.e.,  m  and 
a;  thus,  a  fixed  sampling  routine  with  a  minimal  number  of  points 
could  be  used. 
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Suppose  next  that  X  is  uniformly  distributed  between  a  -  b 
and  a  +  b.  In  this  case  the  characteristic  functions  for  X 
and  S  are  respectively 

*x(u>)  =  (expjaco)  (sin  bu>  r  boo)  (41) 

4>s(u)  =  sin  t  aJTw  (42) 

For  this  distribution  m  =  a,  and  it  is  seen  that  the  modulat¬ 
ing  effect  of  the  non-zero  mean  is  identical  to  that  discussed 
before.  The  characteristic  function  for  the  standardized 
variable  is  again  seen  to  be  independent  of  the  distribution 
parameters  of  X.  In  this  case,  however,  the  function  oscillates, 
and  might  require  more  samples  than  the  Gaussian  case. 

In  both  examples  it  is  noted  that  one  of  the  distribution 
parameters  affects  only  the  displacement  of  the  function, 
and  that  the  other  affects  only  its  spread.  Furthermore,  both 
functions  have  even  symmetry  about  the  mean  value,  with  the 
consequence  that  the  characteristic  function  for  the  standardized 
variables  has  an  imaginary  part  that  is  identically  zero. 

For  the  noise  prediction  model,  the  variable  is  non¬ 
negative,  and  in  general  the  density  function  will  not  be 
evenly  symmetric  about  the  mean.  Thus  the  imaginary  part  of 
the  characteristic  function  for  the  standardized  variable  will 
not  be  identically  zero.  Furthermore,  if  the  ratio  of  o/m  is 
high,  the  density  will  be  high  near  small  positive  values; 
displacing  the  curve  by  the  amount  of  the  mean  value  will 
have  the  undesirable  effect  of  displacing  the  highest  density 
region  from  zero. 
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One  fact  that  can  be  exploited  is  that 

f(x)  =  0,  x  <  0 

Now  f(x)  =  —  /°°  da)  R(co)  cos  ux  +  —  f°°  I  (cc)  sin  cox  (^3) 

11  0  71  0 

where  R(co)  is  the  real  part  of  $  (oi) 

I  (co)  is  the  imaginary  part  of  <t>  (cc) 

For  x  >  0 

f(-x)  =  -■  /  dco  R(co)  cos  uix  —  ~  f  dco  I  (oi)  sin  cox  =  0  (44) 

0  17  0 

Therefore  f(x)  =  -  /“  dco  R(co)  cos  cox,  x  >  0  (45) 

17  0 

_  2  /°°  dco  I  (co)  sin  cox,  x  ^  0  (46) 

ir  0 

That  is,  f(x)  can  be  expressed  as  a  single  integral  for  x  ^  0, 
and  either  (45)  or  (46)  can  be  used. 

This  result  does  not  drastically  reduce  the  computation 
load  since  to  obtain  either  R(co)  or  I(co),  both  real  and  imaginary 
parts  of  the  cumulant  function  must  be  computed;  that  is 
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R  (to) 

=  exp 

S  (id) 

cos  J(to) 

(47) 

I  (to) 

=  exp 

S(to) 

sin  J  (to) 

(48) 

where  S (to)  is  the  real  part  of  the  cumulant  function 

J(io)  is  the  imaginary  part  of  the  cumulant  function 
Substituting  (47)  in  (45)  and  (48)  in  (46)  yields 

2  oo 

f(x)  =  —  /  dto  exp  S(to)  cos  J(io)  cos  tox,  x  >  0  (49) 

11  0 


2  °° 

=  —  /  do)  exp  S  (to)  sin  J  (to)  sin  tox,  x  ^  0  (50) 

0 


Although  (49)  and  (50)  are  equivalent  to  (45)  and  (46)  respective 
the  method  of  their  calculation  could  be  significantly  different. 

Aside  from  the  exponential  spacing  of  the  sample  points, 
the  approach  taken  in  Reference  1  for  the  evaluation  of  (43) 
has  considerable  merit.  In  that  approach,  the  components  of 
the  characteristic  function  are  approximated  by  linear  functions 
between  data  points.  The  resulting  integrals  are  evaluated 
analytically.  This  allows  the  sampling  regime  to  be  dictated 
solely  by  the  components  of  the  characteristic  function.  This 
is  very  desirable  because  the  trigonometric  functions  in  the 
integrands  of  (43)  vary  rapidly  with  to  with  large  values  of  x. 

The  same  approach  can  be  considered  for  the  calculation  of  (49) 
or  (50),  a  practice  that  would  be  desirable  only  if  S(u>)  and 
J(to)  have  a  smaller  "bandwidth"  than  R(to)  or  I ( oj ) .  If  S(to)  and 
J(to)  are  represented  as  linear  functions  between  data  points, 
the  integrals  can  be  evaluated  analytically. 
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An  alternative  that  can  be  considered  is  to  represent  the 
components  as  n—  order  polynomials  in  the  interval  between 
n  +  1  data  points.  With  either  (45)  or  (46)  the  integrals 
would  be  easy  to  evaluate.  The  same  approach  with  (49)  and 
(50)  would  probably  be  intractable  if  n  is  two  or  more. 

At  this  point  the  development  of  (45)  and  (46)  for  n  =  2 

seems  an  attractive  possibility. 

All  of  the  results  of  interest  can  be  restated  for  a 
normalized  variable.  If  it  were  possible  to  determine  in 
advance  that  f(x)  =  0,  x  <  q  >  0,  and  if  q  is  not  small 

compared  to  a,  then  it  would  be  worthwhile  to  compute  the 

characteristic  function  for  the  variable 

T  =  (X  -  q)  *  a  .  (51) 

However,  if  q  cannot  be  determined,  or  if  it  is  small  compared 
to  a,  then  the  reasonable  recourse  is  to  calculate  the  functions 
for 

U  =  X  *  a  (52) 

which  are  ^(w)  =  4>x(io  *  a)  (53) 

and  ^(w)  =  <4x(w  *  a)  (54) 
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The  density  function  for  the  normalized  variable  U  is 

f„(x)  =  —  /°°  du>  R(— )  cos  cox,  x  >  0  (55) 

U  7T  q  O  —  ' 

/°°  du>  1  (— )  sin  x,  x  ^  0  (56) 

0  ° 

From  this  function  various  functions  and  moments  pertaining  to 
X  can  be  calculated: 


_  i 

f X (x)  =  0  fU(x  *  (57) 

XrO 

Fx(x)  =  /  dy  f^y)  (58) 

mn  =  on  /°°  dx  fy(y)xn  (59) 

Methods  are  also  available  for  calculating  Fx(x)  directly  from 
the  components  of  the  characteristic  function.* 

Calculations  for  the  normalized  variable  U  require  a  value 
for  a.  Means  for  its  calculation  or  estimation  are  set  forth 
below,  as  well  as  those  for  the  mean  value. 


♦Albert  H.  Nuttall,  "Alternate  Forms  and  Computational  Considera 
tions  for  Numerical  Evaluation  of  Cumulative  Probability  Distri 
butions  Directly  from  Characteristic  Functions,"  NUSC  Report  No 
NL-3012,  August  1970,  Sec.  2.2. 
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The  derivations  are  based  on  certain  results  given  in 
BBN  Report  No.  3390;  i.e.,  (2-21),  (2-19),  (2-20)  and  (2-22). 

These  results  are  restated  below  for  the  first-order  case  in 
the  same  order.  The  cumulant  function  for  X  is 

m  n 

♦x<“>  =  (60) 

1=1  j=l  1 J 

The  remaining  equations  are 

=  pkij[Hij(u))  -  H..(0)] 

Hi  ■  (w)  =  /  dg  /  dqf  (q)<I>s  (wz) 

R.  1  3 

l 

=  /  dD  /  dr | J (r , D) | f  (q)$  (wz)  (63) 

R.  1  3 

l 

where  is  the  area  defined  by  the  intersection  of  the  sector 
with  the  ith  route.  The  power  transmission  factor  is  implicitly 
a  function  of  the  position  coordinates.  The  mean  and  the  variance 
of  X  can  be  determined  from  the  first  and  second  derivatives  of 
(60)  evaluated  at  zero.  This  involves  the  first  and  second  deri¬ 
vatives  of  (61),  which  in  turn  involves  the  first  and  second 
derivatives  of  (62)  or  (63). 


(61) 

(62) 
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To  obtain  results  for  direct  calculations,  differentiate 
(63)  twice  and  evaluate  for  w  =  0. 

hV  (0)  =  <5>g  (0)  /  dD  /  dr  |  J  (r ,  D)  |  f  (q)z(r,D)  (64) 


Hi .  (0)  =  $s  (0)  f  dD  I  dr|j(r,D)|f  (q)z2(r,D)  (65) 

3  R.  1 

1 

Substituting  (64)  and  (65)  into  the  corresponding  deriva¬ 
tives  of  (60)  and  (6l)  gives  expressions  for  i/v(0)  and  i^(0). 
Exploiting  the  basic  properties  of  cumulant  and  characteristic 
functions  gives 


m  n 

mv  =  p  Z  Z  k..m  /  dD  /  dr  |  J  (r ,  D)  |  f  (q)z(r,D)  (66) 

x  i=i  j=i  Ri 


m  n 

a2  =  p  Z  Z  k .  .  (nu  +  a*  )  S  dD  /  dr  j  J  (r  ,D)  j  f  (q)  z2  (r  ,D) 

X  i=l  j=l  13  Sj  Sj  R  Ui 

Ki 

(67) 


where  m^  is  the  average  of  the  mean-squared  pressure 
Sj  of  sources  of  type  j 

a2  is  the  variance  of  the  mean-squared  pressure 
of  sources  of  type  j. 
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Each  term  of  (66)  or  (67)  can  be  regarded  as  an  extension 
of  Campbell’s  theorem*,  which  pertains  to  moments  of  elemental 
shot  noise. 

To  obtain  results  for  the  approximations,  differentiate 
(62)  twice  and  evaluate  for  00  =  0: 


H± j  (0 )  “  $s  .  f  dq  -T  dqfQ  (q)  z  (r , D)  (68) 

3  R.  i 

1 

(°)  =  (0)  /  dg  f  dqf  (q) z2  (r,D)  (69) 

j  wi 

If  the  power  transmission  factor  does  not  vary  appreciably  over 
R^,  then 

z(r,D)  =  zi  (70) 

where  is  the  power  transmission  factor  from  the  reference 
point  in  to  the  sensor. 

If  (70)  is  substituted  in  (68)  and  (69),  the  integrals 
can  be  readily  evaluated,  and  the  results  are 


*A.  Papoulis,  Probability,  Random  Variables,  and  Stochastic 


Processes ,  P.  359. 
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Hij(0)  *  $s.(0)^izi  (71) 

H- j(°)  *  (O)^z^  (72) 

where  is  the  length  of  the  route  in  R^* 

Substituting  (71)  and  (72)  into  the  corresponding  derivatives 
of  (60)  and  (61)  gives  approximate  expressions  for  \p'  (0)  and 
*J(0).  Exploiting  the  basic  properties  of  cumulant  and  charac¬ 
teristic  functions  gives 


m 


X 


m  n 
p  Z  l 
i=l  j=l 


(73) 


m  n 

l  E  k .  .  (mi 
i=l  j=l  13  bj 


+  a 


s.]kzl 


(74 ) 


The  quantity  Is  the  average  number  of  ships  of  type  j  in 

V 

The  problem  of  determining  the  points  at  which  to  calculate 
the  characteristic  function  would  be  simple  if  the  function  were 
known  in  advance.  If  that  were  the  case  of  course,  there  would 
be  no  need  to  calculate  it.  However,  it  will  be  possible  to 
obtain  an  analytic  approximation  if  (70)  holds,  and  if  the 
source  characteristic  function  is  analytic.  Substituting  (70) 
in  (62)  and  evaluating  the  integrals  gives 
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,  ^  _ 


Hi j  (aj)  =  Uzi)  ( 75 ) 

j 

and  substituting  this  result  in  ( 6 1 )  gives 

^(w)  =  [4>s  (o)Zi)  -  1]  (76) 

Finally,  substituting  this  result  in  ( 6 0 )  gives  the  cumulanl 
function  for  X.  However,  it  may  not  be  necessary  to  do  a 
complete  evaluation.  It  appears  from  (76)  that  the  critical 
routes  are  those  with  the  larger  power  transmission  factors  z.. 
Also,  it  might  be  that  the  source  characteristic  function  for 
one  class  of  ships  has  more  detail  than  the  others  and  necessi¬ 
tate  more  closely-spaced  points  uk.  The  result  should  be  a  good 
approximation  for  sector  noise,  but  it  will  not  be  as  accurate 
for  beam  noise. 

3.2  An  Inverse  Fourier  Transform  Algorithm 

The  purpose  of  this  section  is  to  develop  the  equations 
for  an  approximation  technique  to  evaluate  the  special  form 
of  the  inverse  Fourier  transform  given  by  (45)  of  Section  3.1: 

f(x)  =  -  f+°°  R(w)  coswx  dw,  (77) 

11  0 

for  x  >  0;  where  f(x)  is  a  probability  density  function  such 
that  f(x)  =  0,  for  x  <  0.  R(w)  is  the  real  part  of  the  charac¬ 
teristic  function  of  f(x). 
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The  first  assumption  is  that  R(w)  can  be  approximated  by 
a  quadratic  in  w  over  a  three-sample-point  interval,  that  is 

R(w)  =  aco2+bw  +  c  (78) 

for  the  general  interval  £  w  £  t0n+2*  The  constants  a,  b  and 
c  can  be  determined  from  the  set  of  equations 

R(n)  =  a  a>n2  +  b  un  +  c 

R(n+1)  =  a  ojn+12  +  b  un+1  +  c  (79) 

R(n+2)  =  a  u)n+22  +  b  u)n+2  +  c 

where  R(n)  is  the  value  of  the  real  part  of  the  characteristic 

function  at  to  =  to  ,  etc. 

n 


R(«) 


FIGURE  5  Approximating  Polynomial 
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If  it  is  also  assumed  that  the  samples  are  equally  spaced 
in  frequency  w,  then 

u>n+2  =  wn+1  +  Aw  =  a)  +  2  Aw  ( 8  0 ) 

where  Aw  is  the  interval  between  samples  (See  Figure  5.) 

The  solutions  for  the  constants  a,  b  and  c  in  terms  of 
Aw,  R(n),  R(n+1)  and  R(n+2)  are  given  by 


a  =  — -  [R(n)  -  2R(n+l)  +R(n+2)] 

2Aw2 


(81) 


b  =  -=L 


2Aw' 


|  2wn  [R(n)  - 


2R(n+l)  +  R(n+2 ) ] 


+  Aw  C3R(n)  -  4R(n+l)  +  R(n+2)]  >  (82) 


2Aw' 


jlKn)  (wn 


+  2  Aw)  (wn  +  Aw  ) 


2<»>n  R(n+1)  (u^  +  2  Aw  ) 


+  w  R(n+2)(w  +  Aw  ) 
n  n 


(83) 


Consider  the  inverse  Fourier  integral  of  Equation  (77)  eval¬ 
uated  over  the  general  interval  wfi  to  wn+2,  for  some  value  of 
x  >  0: 


i,\  ’ 1 1  rir ' ' r  ~ 


. _ 
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Afn(x) 


/ 


.“n+2 

R(<*>)  cos  ux  dw 


(i) 


n 


Using  the  quadratic  approximation  for  R(o>)  given  in  Equation 
(78),  we  have 


0) 


“.w  5 


/ 


u>_ 


n+2 


[ao>2  +  b  a)  +  c]  cosidx  dw  (84) 


n 


where  the  constants  a,  b  and  c  are  given  in  Equations  (81) 
through  (83).  This  integral  can  be  evaluated  exactly: 


wn(x) ; 


03=0) 


2aM-|-b  cos  «  +  sin  „x 


n+2 


03=03 


n 


(85) 


Evaluating  Equation  (85)  at  the  upper  and  lower  limits  and 
substituting  the  expressions  for  a,  b  and  c  from  Equations  (8l) 
through  (83)  gives 


if  (x)  =  B.(n)-Wnti)U.»(nrt),  oos  „ 
n  2  iu  x2  n+2 


1  .  .  . .  ■  .iV-a 
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FIGURE  6  Sampling  Scheme 


With  these  assumptions,  the  approximation  for  f(x)  can  be  simplified 
after  considerable  algebraic  manipulation 

f(x)  ^  3R(0)  -  ^R(l)  +  R(2) 
irAw  x  2 


- -  NE  |  [R(n-2)  -  ^R(n-l)  +  6R(n)  -  ^R(n+1)  +  R(n+2)1 

iU)  x  2  n=2  l 


cosu  X 
n 


(n  goes  by 
steps  of  2) 
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2R(n-2)  -  ^R(n-l)  +  HR(n-H)  -  2R(n+2) 


sin  w 


+  R(N-2)  -  4R(N-1)  +  3R(N) 
irAw  x  2 


cos  wXTx 
N 


2  2 


+  2R(N)x  -^2[R(N-2)  -2R(N-1)  +R(N)] 


sin  u).x  x 
N 


for  x  >  0 


A  recursion  relationship  can  be  derived  for  sin  w  x  and  cos  w  x 

n  n 

because  of  the  equal  spacing  of  samples  in  the  frequency  domain. 

Note  that  to  =  nAw,  where  Aw  is  the  interval  between  samples, 
n 

Therefore, . 


sin  w  x  =  sin  nAwx 
n 


s  2  sin  (n-l)Awx  cos  Aiox  T  sin  (n-2)Awx  (88) 


ill 
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The  recursion  relation  becomes 

sin  a>nx  =  [2  cos  Awx  ]  sin  to  ^  x  -  sin  u>n  ^  x  (89) 

The  first  three  values  are 

sin  w  x  =  0 
o 

sin  w  x  =  sin  Aw  x 
1 

sin  w^x  =  [2  cos  Aux]  sinAwx 
The  recursion  relationship  for  cos  wnx  is 

cos  o>nx  =  [2  cos  Aux]  cos  wn_ix  “  cos  u)  n_2x  (90) 

The  first  three  values  are 


cos  wqx 


cos  w  x 

1 

COS  0)  X 

2 


1 

cosAwx 

2[cos  Awx  ]2  -1 


The  numerical  approximation  formula  for  f(x)  given  by 
Equation  (87)  is  valid  for  any  x  >  0.  For  the  case  where 
x  =  0,  we  must  use  a  different  formula.  From  Equation  (77) 
we  have 


f (0)  =  - 


R(u>)  dm 


(91) 
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Using  the  quadratic  approximation  for  R(w)  given  by  Equations 
(78)  and  (8l)  in  Equation  (91)  >  f(0)  is  approximately  equal  to 


f(0)  =  ^  [R(N)  -R(0)] 


N-2 

r  [R(n)  +  2R(n+l)] 
n=0 


(92) 
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4.0  PROGRAM  INPUTS 

The  prediction  of  the  beam  noise  density  function  requires 
acoustic  data  pertaining  to  the  sensor,  its  acoustic  basin,  and 
surface  ships.  Section  4.1  discusses  a  beam  pattern  function 
for  a  shaded  line  array  of  hydrophones.  Section  4.2  gives  route 
geometry  data  for  the  Northeast  Pacific  Ocean,  and  Section  4.3 
presents  required  source  characteristic  functions  for  classes 
of  surface  ships. 

4.1  Beam  Pattern  Function 


In  the  calculations  described  in  Reference  1,  only  the 
noise  in  a  narrow  sector  due  to  merchant  shipping  was  considered. 
For  the  present  effort,  the  entire  beam  pattern  is  to  be  incor¬ 
porated  into  the  calculations.  The  approach  given  in  the  sequel 
is  to  use  the  normalized  beam  pattern  as  a  weighting  coefficient 
to  be  applied  for  a  given  angle  in  the  calculations  <J>Y. 


Engineers  of  the  Naval  Research  Laboratory  supplied  a 
beam  pattern  approximating  Hanning  shading: 


B(  0 ) 


2 


N  sin  tt g 
2  it  3  ( 1— 3 2 ) 


(93) 


where 

3=  (sin  9  -  sin  0S) 
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and 


0  =  beam  steered  direction 

s 

0  =  angle  from  broadside 
L  =  length  of  the  array 
N  =  number  of  phones 
X  =  c/f 

c  =  speed  of  sound 
f  =  frequency 


Figure  7  shows  the  reference  for  the  array  angles. 


FIGURE  7  Array  Angles 
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Let  f  equal  the  design  frequency  and  cq  equal  the  design 
value  of  the  speed  of  sound.  The  phone  spacing  d  is  given  by 

X.  c 


length  of  the  array  L  is  given  by 


L  =  (N-l)d  =  (N-l )  2Y~ 

o 


Define  the  constant  a  =  — .  Therefore, 

X 


_  (M-l)  V 
a  2  cf 

o 


(96) 


Table  2  gives  the  value  of  a  for  four  selected  frequencies 
A  design  speed  of  sound  cq  =  5 ,000  ft. /sec.  and  an  actual  speed 
of  sound  of  c  =  4,860  ft. /sec.  were  used  in  the  computation  of 

a. 


Table  2:  Values 

of  a  for 

Four  Frequencies 

f(Hz) 

fo(Hz) 

N 

a 

10 

20 

62 

15.68930041 

40 

60 

58 

19.5473251 

50 

60 

58 

24.43415638 

300 

320 

64 

30.38194444 

>  iirti  im  ii  i  ■  r  ir  n  r  i  aiill 
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The  normalized  beam  pattern  BP  is  given  by 

BP  °[rrvi!e>)]  ‘  <97> 

where  3  =  a  (sin  6  -  sin  0  ). 

s 

Note  that  there  are  three  points  (3=0  and  3  =  +1 )  at  which 
the  function  is  an  indeterminate  form.  Using  L' Hospital's  rule 
BP  can  be  evaluated  at  these  points 

BP  =  1.0  at  3  =  0  (98) 


and 


BP  =  0.25  at  3  =  +1  (99) 

Figure  8  shows  the  normalized  beam  pattern,  converted  to 
dB  for  a  frequency  of  10  Hz  and  a  steering  angle  0g  of  -9 
degrees  (351°). 

Figure  9  shows  the  normalized  beam  pattern,  converted  to 

dB  for  a  frequency  of  50  Hz  and  0  =  -9  degrees.  Note  that 

s 

as  the  frequency  increases,  the  main  lobe  becomes  narrower  and 
there  are  more  side  lobes. 

In  the  computations  of  the  characteristic  function  for  Y, 
there  is  an  integral  over  the  angle  D.  The  step  size  for  this 
intergral  will  have  to  be  adjusted  such  that  the  time  for  com¬ 
putation  on  the  digital  computer  is  not  excessive.  One  approach 
is  to  use  a  finer  step  size  for  angles  in  the  main  lobe,  and  a 
coarser  step  size  for  angles  in  the  side  lobes.  However,  for 
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the  side  lobes  an  approximation  must  be  used  that  is  compatible 
with  coarser  angle  step  sizes.  For  the  side  lobes,  the  norma- 
lized  beam  pattern  will  be  modelled  as  one-half  the  envelope 
value,  that  is 


BP  =  1/2 


[  ttB  U-B2)  ] 


(100) 


Note  that  the  main  lobe  occurs  for  -  2  ^  3  <  +  2  and  the 
side  lobes  occur  for  values  of  3  outside  this  range. 

In  Figures  8  and  9  the  one-half  envelope  approximation 
is  shown  as  the  solid  line  through  the  side  lobes. 


4.2  Route  Geometry  Data  for  Northeast  Pacific  Ocean 


This  section  sets  forth  route  data  for  the  Northeast 
Pacific  Ocean. 

Eight  routes  were  selected  from  the  seventy-five  canon¬ 
ical  routes  described  in  Reference  5.  However,  in  some 
instances,  portions  of  the  canonical  routes  are  eliminated 
because  it  is  felt  that  they  contribute  little  to  the  noise 
statistics.  For  example,  it  is  felt  that  the  Aleutian  chain 
forms  an  effective  boundary  for  the  Northern  Pacific;  therefore, 
the  canonical  route  from  Honolulu  to  Nome  is  terminated  at  the 
Unimak  Pass.  Further,  due  to  internal  constraints  of  the  pro¬ 
gram  module  GEOM* ,  two  of  these  eight  routes  must  be  counted 
as  double  routes.  Therefore,  the  example  contains,  in  effect, 
ten  routes. 


^Described  in  Section  VII  of  Reference  6. 
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Each  route  consists  of  one  or  more  segments,  with  the 
following  items  needed  for  each  segment: 

1.  latitude  and  longitude  of  the  end-points 

2.  width  at  the  end-points 

3.  maximum  width 

4.  ship  densities  for  each  class  of  vessel 

Table  3  contains  Items  1,  2  and  3  for  each  of  the  ten  routes 
considered  in  this  case.  Route  and  segment  numbers  are  de¬ 
scriptors  used  by  the  program  GEOM.  In  addition,  geographic 
designations  are  given  for  each  segment  end-point.  Figure  10 
is  a  map  of  the  Pacific  Ocean,  indicating  (approximately)  the 
positions  of  the  routes. 

4.3  Source  Characteristic  Functions 

It  is  shown  in  Section  2.1  that  the  calculation  of  the 
characteristic  function  of  Y  requires  a  source  characteristic 
function  $Sj(  )  for  each  type  of  ship.  In  particular,  each 
of  these  characteristic  functions  pertains  to  the  same  frequency 
band  as  the  ambient  noise  prediction. 

The  required  source  characteristic  functions  were  derived 
from  statistical  information  produced  by  Heine  and  Gray  [31 
pertaining  to  principal  sources  of  noise  radiated  by  merchant 
ships.  This  information  was  employed  to  derive  discontinuous 
probability  density  functions  for  the  level  of  each  type  of 


source. 
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SEGMENT  ORIGIN 


ROUTE 

NUMBER 


NUMBER 

OF 

SEGMENTS 


SEGMENT 


LATITUDE 

LONGITUDE 

45°43 

•N 

1 42°40 

'E 

54o20 

•N 

1 65ol*0 

'W 

35°27 

•N 

1 39°39 

■E 

43“  0 

•N 

150°  0 

'E 

2*3°  0 

*N 

170°  0 

'W 

43°  0 

'N 

160°  0 

'W 

43“  0 

'N 

160°  0 

'W 

21  °40 

'N 

157°30 

'W 

21  °40 

'N 

1 57°30 

’W 

2  1  °40 

'N 

157“30 

'W 

6°55 

'N 

81  °25 

'W 

37°48 

'N 

122“24 

'W 

54°20 

'N 

165°  0 

•w 

SEGMENT 

WIDTH 

(Degrees 


.25 


.125 


GEOGRAPHIC 

DESIGNATION 


La  Perouse 
Stra  i  t 
(East  End) 


Un i mak 
Pass 


Yokohoma 


Mid-Ocean 


M id-Ocean 


M id-Ocean 


Mid-Ocean 


Honolulu 


48  30' 
I24“40' 

40° 1 8’ 
I24“30' 


Honol u 1 u 

37°35'N 

1 22°40'W 

Isla 

Coiba 

21 °30 1 N 
157“35'W 

San 

Francisco 

17°32'S 

149°34’W 

Unimak 

Pass 

48°30 1 N 
124°40'W 

Strait  of 
Juan 

De  Fuca 

40o18'N 

124o30'W 

Cape 

Mendoc i no 

37°48'N 

122°24'W 

SEGMENT  TERMINUS  SEGMENT 

SEGMENT  MAXIMUM 

LATITUDE  WIDTH  GEOGRAPHIC  WIDTH 

LONGITUDE  (Degrees)  DESIGNATION  (Degrees) 


Unimak 

Pass 


Strait  of 
Ju 

De  Fuca 


Mid-Ocean 


Mid-Ocean 


M i d-Ocean 


Strait  of 
Juan 
De  Fuca 


San 

Francisco 


Unimak 

Pass 


Strai t  of 
Juan 
De  Fuca 


San 

Francisco 


Honolulu 


Papeete 


Strait  of 
Juan 
De  Fuca 


Cape 

Medoc i no 


San 

Francisco 


54° 

20 

N 

165° 

0 

W 

48° 

30 

N 

124° 

40 

W 

43° 

0 

N 

150° 

0 

E 

43“ 

0 

N 

170° 

0 

W 

43° 

0 

N 

160° 

0 

U 

48° 

30 

N 

124° 

40 

U 

37° 

48 

N 

122° 

24 

W 

54“ 

20 

N 

165“ 

0 

'W 

TABLE  3 


ROUTE/ SEGMENT  GEOMETRIC  DATA 
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The  procedure  for  obtaining  the  source  characteristic 
function  from  the  discontinuous  density  functions  is  out¬ 
lined  below. 


1.  Select  a  characteristic  function  with  at  least 
two  parameters  that  applies  to  a  non-negative 
random  variable. 


Calculate  the  low-order  moments  of  the  source 
intensity  of  each  component.  These  can  be 
readily  obtained  from  the  discontinuous 
density  functions  fy  (  )  for  the  component 
source  levels.  For  example 5  the  mean  value 
of  the  source  intensity  of  the  i^  component 
is 


mi  =  /  du  10u:1°  fU;L(u)  (101) 

3.  Determine  the  low-order  moments  of  the  source 
intensity  from  the  low-order  moments  of  the 
component  source  intensities. 

4.  Determine  the  values  of  the  parameters  of  the 
characteristic  function  so  that  its  moments 
match  the  calculated  moments. 

The  advantage  of  this  procedure  is  that  it  avoids  the  excessive 
oscillations  of  the  characteristic  functions  that  would  be  in¬ 
troduced  by  the  direct  utilization  of  the  discontinuous  density 
functions. 
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The  remainder  of  this  section  (1)  summarizes  the  informa¬ 
tion  developed  by  Heine  and  Gray,  (2)  presents  plots  of  the 
probability  density  functions  derived  from  them,  and  (3)  gives 
the  low-order  moments  of  the  source  intensity  derived  from  the 
source  level  density  functions. 

Heine  and  Gray  developed  histograms  for  the  source  level 
and  frequency  for  four  categories  of  radiated  noise.  The 
spectral  components  considered  were  a  broadband  component  due 
to  propeller  cavitation;  lines  at  the  propeller  blade  rate  and 
harmonics;  lines  at  the  diesel  firing  rate  and  harmonics;  and 
a  line  at  the  electric  plant  frequency.  The  histograms  pertain 
to  "Long  Range  Effective  Source  Level"  (LRESL),  defined  as  [31 
the  power  radiated  by  an  omnidirectional  source  (at  a  depth  of 
20  feet)  between  depression  angles  zero  and  fifteen  degrees 
(including  that  from  surface  interference) . 

4.3.1  Broadband  Cavitation 


Heine  and  Gray  found  it  useful  to  divide  the  merchant  fleet 
into  three  classes  of  vessel-fishing,  small  merchant,  and  large 
merchant.  In  addition,  each  of  these  was  divided  into  eight 
acoustic  classes.  Histograms  showing  the  distribution  of  acoustic 
class  for  each  class  of  vessel  are  shown  in  Figures  11,  12  and  13, 
while  their  broadband  spectra  are  shown  in  Figure  14. 

Assuming  statistical  independence  between  the  acoustic 
class  and  broadband  spectrum,  composite  histograms  are  formed 
for  the  source  level  at  the  frequencies  10,  40,  50  and  300  Hz 
(in  a  1-Hz  band),  for  the  three  classes  of  vessel.  Finally, 
by  considering  10  dB  increments  and  normalizing  the  area  to 
unity,  an  approximate  probability  density  function  is  generated 
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FIGURE  11  Distribution  of  Broadband  Classes,  Fishing  Vessels 
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FIGURE  12  Distribution  of  Broadband  Classes,  Merchant  Vessel 
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FIGURE  13  Distribution  of  Broadband  Classes,  Merchant  Vessel 

Greater  than  700  Feet  LOA 


Broadband  Acoustic  Class  Spectra 
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from  each  composite  histogram.  These  density  functions  are 
shown  in  Figures  15  through  26 

4.3.2  Propeller  Blade  Rate  and  Harmonics 

Heine  and  Gray  I3J  developed  histograms  for  the  source 
level  and  the  fundamental  frequency  for  the  three  classes  of 
vessel.  Gray  [6]  modified  the  source  level  distribution  and 
developed  a  [7]  model  for  the  effects  of  harmonics. 

Assuming  statistical  independence  between  source  level 
and  frequency,  probability  density  functions  were  generated 
for  source  level  at  the  frequencies  10,  40,  50  and  300  Hz  (in 
a  .2-Hz  band)  for  the  three  classes  of  vessel.  These  distribu¬ 
tions  include  the  effects  of  the  first  nine  harmonics.  Note 
that  fishing  vessels  contribute  nothing  at  10  Hz,  while  small 
and  large  merchant  contribute  nothing  at  300  Hz.  These  mono¬ 
pole  source  level  density  functions  are  shown  in  Figures  27 
through  35-  The  figures  give  values  for  P(-°°),  which  is  the 
probability  that  there  is  no  line  in  the  specified  1-Hz  band. 
The  procedure  employed  is  described  in  Section  4.2  of  Refer¬ 
ence  1. 

4.3.3  Electric  Plant  Frequency 

Heine  and  Gray  [3]  developed  histograms  for  the  source 
level  and  frequency  for  the  three  classes  of  vessel.  Assum¬ 
ing  statistical  independence  again,  a  probability  density 
function  was  generated  for  the  source  level  at  5C  Hz  for  each 
class  of  vessel.  Since  the  frequencies  that  are  of  interest 
are  10,  40,  50  and  300  Hz,  the  line  at  60  Hz  does  not  contribut 
any  power.  The  probability  density  functions  are  shown  in 
Figures  36,  37  and  38. 
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FIGURE  18  Broadband  Source  Level  Distribution  for  Fishing  Vessels  at  300  Hz. 


Bolt  Beranek  and  Newman  Inc 


Report  No.  3653 


65 


FIGURE  20  Broadba 


Broadband  Source  Level  Distribution  for  Small  Merchant  Ships  at  50  Hz- 


FIGURE  22  Broadband 
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FIGURE  29  Blade  Rate  Source  Level  Density;  Fishing  Vessel,  300  Hz. 
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FIGURE  32  Blade  Rate  Source  Level  Density;  Small  Merchant  Ship,  50  Hz. 


Report  No.  3653 


i 


— l - 

o 

Cn! 

Q  Ajjijqeqojj 


1A 

d 

\D 


o 


FIGURE  33  Blade  Rate  Source  Level  Density;  Large  Merchant  Ship,  10  Hz. 
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FIGURE  35  Blade  Rate  Source  Level  Density;  Large  Merchant  Ship,  50  Hz. 
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FIGURE  36  Narrowband  Electric  Plant  Source  Level  Density;  Fishing  Vessel,  50  Hz. 
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Revised  information  concerning  the  diesel  firing  rate 
and  its  harmonics  are  not  yet  available. 

4.3.4  Means  and  Variances 


The  parameters  required  to  establish  the  source  charac¬ 
teristic  functions  are  given  in  Tables  4  through  6,  which 
gives  means  and  variances  for  rms  pressure  in  a  0.2  Hz  band 
for  monopole  sources. 


Table  4:  Mean  and  Variance  for  Broadband  Source 

Large  Merchant 

Small  Merchant 

Fishing 

Frequency  (Hz) 

m  (yP)2 

1.81 06x1 0 1 5 

6.6383xl014 

4.4093xl013 

10 

a2  (yP)4 

5.6936X1031 

9. 3840x1 03° 

6.7276xl027 

10 

m 

1 . 1 320xl016 

4.1505xl0ls 

2.5329xl014 

40 

a2 

9.1 073x1 0 3  2 

3.6672xl032 

2. 6288x1 02 9 

40 

m 

7 . 241 6x1 01 5 

2.6552xl0ls 

1.6204x10s 4 

50 

a2 

9.1 073x1 032 

1.5010x10s2 

1.0759X1029 

50 

m 

8. 8897x1 09 

5. 541 7x1 09 

3. 4314xl09 

300 

a2 

2. 401 1 xl 02  8 

1. 8089x1 02 8 

8.4460xl027 

300 

*Note:  Monopole  source  rms  pressure  in  a  0.2  Hz  band. 
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Table  5:  Mean 

and  Variance  for 

Blade  Rate  Source 

Large  Merchant 

Small  Merchant 

Fishing 

Frequency  (Hz) 

m  (yP)2 

2.996xl018 

3.643xl017 

0.0 

10 

o2  (yP)4 

3.028x10 

1.637xl038 

0.0 

10 

m 

2.275x10 

3.847X1011* 

2. 800x1 01 5 

40 

o2 

2.223x10 

1 . 61 2x1 032 

1 . 528xl035 

40 

m 

1.633x10 

1.035xl014 

7.850xl014 

50 

o2 

7.061x10 

1.456xl031 

1.872xl034 

50 

m 

0.0 

0.0 

9. 800x1 08 

300 

o2 

0.0 

0.0 

9.783xl025 

300 

♦Note:  Monopole  source  rms  pressure  in  a  0.2  Hz  band. 


Table  6: 

Mean  and  Variance,  Electric  Plant  Line  ac  ou  nz 

Large  Merchant 

Small  Merchant 

Mean  (yP ) 2 

1 . 3671 xl 0 1 5 

3.141 9x1 0 1 4 

4.2394xl013 

Variance 

(yP)“ 

7. 0008x1 031 

7.8622X1029 

7. 7548x1 02 8 
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5.0  PREDICTION  OF  NARROWBAND  COMPONENTS 

The  noise  resulting  from  the  operation  of  surface  ships 
is  not  unlike  that  resulting  from  the  operation  of  threat  sub¬ 
marines.  In  particular,  the  presence  of  narrowband  components 
(lines)  in  the  radiated  noise  signatures  of  surface  ships  makes 
the  acoustic  detection  of  threat  submarines  more  difficult. 
Predicting  the  effects  of  lines  on  the  performance  of  surveill¬ 
ance  systems  requires  either  statistics  or  statistical  measures. 
One  measure  of  interference  is  the  number  of  lines  (1)  whose 
center  frequencies  fall  within  a  specified  frequency  interval 
Af,  and  (2)  whose'  intensities  at  the  point  of  observation  are 
greater  than  a  specified  threshold  value  t.  If  the  threshold 
value  t  is  appropriately  selected,  the  measure  is  observable 
by  means  of  narrowband  analysis  equipment. 

The  analysis  which  follows  pertains  to  interference 
generated  by  ships  in  a  specified  sector  of  azimuth.  Exten¬ 
sion  to  array  beam  interference  is  discussed  later. 

The  total  number  of  lines  satisfying  the  enumerated 
criteria  is  a  random  variable  that  can  be  expressed  as  the 
sum  of  contributions  from  individual  ships: 


N  = 


A ,  , 

m  n  ij 
III 
i=l  j=l  k=l 


(102) 
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where  m  is  the  number  of  route  envelopes  intersecting 

the  sector 

n  is  the  number  of  ship  types 

A.,  is  the  number  of  ships  of  type  j  on  route  i, 

0  a  random  variable 

N.  ..  is  the  number  of  lines  satisfying  the  criteria 
radiated  by  ship  k  of  type  j  on  route  i. 

Each  ship  has  several  mechanisms  that  radiate  narrowband  funda¬ 
mentals  and  harmonics;  their  effect  is  tallied  by 


p  » 

Z  Z 
1=1  h=0 


Uijki?h 


(103) 


is  the  number  of  types  of  sources 

is  a  random  variable  pertaining  to  har¬ 
monic  h  of  source  /  on  ship  k  of  type  j 
on  route  i;  its  value  is  one  if  the 
enumerated  criteria  are  satisfied,  and 
zero  if  not. 

The  random  variable  j k:^ki  can  be  expressed  as  a  function 
of  the  characteristics  or  the  source  and  of  the  medium: 


uijkih=  0(pijk/h)  u  IPijk/h z  (Gijk>  QiJk 


where  is  the  center  frequency  of  the 

^  **  designated  by  the  subscripts 

P.  ..  n,  is  the  source  .intensity  of  the 
^  *  designated  by  the  subscripts 


-  t]  (104) 

line 

line 


where  p 

Ui j  k/h 
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G. ..  is  the  longitudinal  (along  the  route) 

1*'  coordinate  of  the  position  of  ship  k 
of  type  j  on  route  i;  it  is  a  random 
variable  that  is  statistically  independent 
of  the  longitudinal  coordinate  of  any 
other  ship. 

Q. ..  is  the  transverse  coordinate  of  the 

^  position  of  ship  k  of  type  j  on  route  i; 

it  is  a  random  variable  that  is  statistically 
independent  of  all  other  position  coordinates 


e(f )  =  1,  if  f  is  in  Af 


=  0,  if  not 


This  function  counts  all  of  the  lines  in  the  designated 
frequency  interval. 

u(x)  =  1,  x  >  0 

=  0,  if  not 

This  function  counts  all  lines  whose  intensities  are  above 
the  designated  threshold  at  the  observation  point. 


Various  statistical  measures  of  the  line  count  can  be 
derived  using  (102),  (103)  and  (104).  For  this  initial 
effort,  the  ensemble  mean  or  expected  value  will  be  derived. 
This  measure  is  one  that  has  been  estimated  in  experimental 
work.  The  mean  value  of  the  line  count  is 
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Substituting  (102)in  (105)  and  performing  the  expec¬ 
tation  yields 


m  n  00  a 

mM  =  E  E  E  p.  .  (a)  E  E(N.  ..  )  (106) 

N  i=l  j=l  a=0  1J  k=l  ijk 


where  P-^j  (  )  Is  the  frequency  function  for  A jj  .  Since 
E(Nijk)  does  not  depend  on  the  distribution  parameter  a 
[See  (3)1,  all  terms  of  the  inner  sum  are  equal,  and 


mN  = 


m 

E 

i=l 


n 

E  E(N 

j=l 


ijk; 


E  p  .(a)a 
a=0 


(107) 


The  inner  sum  of  (107)  is  the  ensemble  average  A.^  of  A 
thus. 


m  n 


-M  =  Alj  E(Nijk) 


The  expected  value  of  (103)  yields 

P  00 


**  *  ifi  Jo  E(lW 


(108) 


(109) 


and  the  expected  value  of  (104)  can  be  stated  as 
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OO  00 

/  df  f  (f)e(f )  /  dp  fp  (p)  /  dg  /  dq 

0  Fdh  0  Pj£h  p 

i 


x  fr  (g)f0  (q)u[pz(g,q)  -  t]  (110) 

yi 


where  R^  denotes  a  region  which  includes  the  intersection  of 
the  observation  sector  and  the  envelope  of  route  i. 


It  will  be  assumed  that  the  distribution  of  the 
nates  Gy^  is  uniform;  furthermore,  a  change  from  the 
variables  to  the  observation  variables  simplifies  the 
tion  of  Ry  With  these  changes,  the  result  is 


coordi- 

route 

descrip- 


E(U...  ah)  =  sT1  /  dffp(P)  f  dPfP  ,  f  dD  1  dr I J ( x* , d )  | 
ijkUh  i  Af  o  r 


x  fn  (q)u[pz(r,D)  -  t] 


(ill) 


where  gr1  is  the  arc  length  of  the  route  in  Ry 

Substituting  (111)  into  (109)  and  the  result  into  (108) 

gives 


m  n 

nijj  ®  p  Z  Z  k 


Z  Z  k  Z  Z  /  df  f  ■  (f)  /  dpfp  (p) 
i=l  j=l  0  h=0  Af  Fjih  0  Pj/h 


x  /  dD  /  dr|j(r,D)|fQ  (q)u[pz(r,D)  -  t] 
*1 


(112) 
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where  p  Is  the  earth's  radius 

k..  =  A. .  t  (pgi)  is  the  average  number  of 
J  J ships  of  type  j  on  a  unit  length 
of  route  i. 

Since  z(r,D)  is  zero  outside  of  the  sector,  it  is  clear  that 
the  integration  on  D  can  be  restricted  to  the  sector.  The  rela¬ 
tionship  between  the  route  coordinates  and  the  observation 
coordinates  are  given  in  Section  2.2  of  Reference  1. 

Reference  3  presents  information  on  the  narrowband  source 
characteristics  of  three  classes  of  ships:  fishing,  and  small 
and  large  merchant.  Histograms  are  given  for  the  center  frequency 
and  source  intensity  generated  by  the  propeller,  the  electric 
plant,  and  by  the  propulsion  machinery.  The  histograms  apply 
directly  to  the  fundamental  output,  and  can  be  modified  to 
apply  to  harmonics.  For  the  frequency  histogram,  the  abscissa 
values  are  multiplied  by  (h  +  1),  where  h  is  the  harmonic  number. 
Equations  are  given  that  relate  the  source  level  of  a  harmonic 
to  that  of  the  fundamental.  The  difference  between  the  two  is 
the  amount  by  which  the  abscissa  values  are  changed. 

One  could  consider  making  Doppler  corrections  for  the 
frequency  histograms  according  to  frequency,  speed,  and  direc¬ 
tion.  For  most  applications,  the  correction  will  probably  not 
be  required. 

Equation  (112)  pertains  to  the  average  number  of  lines 
received  at  a  point  that  lie  within  a  designated  frequency  band, 
that  are  above  a  specified  threshold,  and  that  are  radiated  by 
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G 

ships  in  a  designated  azimuth  sector.  The  result  is  readily 
modified  for  the  prediction  of  array  beam  noise.  For  this 
case,  the  argument  of  the  unit  step  function  becomes 

pz (r,D)b [D,E(r ) ]  -  t  (113) 

where  b(D,E)  is  the  power  gain  function  of  the  array- 
beamformer  combination  for  a  plain-wave 
arrival  in  terms  of  the  azimuth  and  ele¬ 
vation  angles. 

The  region  of  integration  includes  the  entire  azimuth 
sector  between  the  points  at  which  the  route  enters  and  leaves 
the  acoustic  basin.  This  accounts  for  the  contributions  of 
the  sidelobes  as  well  as  the  main  lobe. 
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